import numpy as np
import matplotlib.pyplot as plt

# Set up the plot configurations.
plt.rcParams['lines.linewidth'] = 2.0
plt.rcParams['lines.color'] = 'black'
plt.rcParams['legend.frameon'] = True
plt.rcParams['figure.figsize'] = (8, 6)
plt.rcParams['font.family'] = 'serif'
plt.rcParams['legend.fontsize'] = 15
plt.rcParams['font.size'] = 15
plt.rcParams['axes.spines.right'] = False
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.left'] = True
plt.rcParams['axes.spines.bottom'] = True
plt.rcParams['axes.axisbelow'] = True
plt.rcParams['grid.color'] = 'grey'
plt.rcParams['grid.alpha'] = 0.6
plt.rcParams['grid.linestyle'] = '--'

sizex = 7.5
sizey = 5

# Create a figure and axis.
fig, ax = plt.subplots(2, 1,figsize=(sizex,2*sizey))

F=1
E=1
A=1
l=20
n=1

x = np.linspace(0, l, 100)


# Plotting the functions.
ax[0].plot(x, (F/(E*A)+n*l/E)*x-(n/(2*E))*x**2,color='red',linewidth=2.0,label="analytische Lösung")

fem_sol=np.array([
 [0,0],
 [l/2,F*l/(2*E*A)+3*n*l**2/(8*E)],
 [l,F*l/(E*A)+n*l**2 / (2*E)]   
])

#ax[0].scatter(fem_sol[:,0],fem_sol[:,1],color="blue",marker="o",s=100,label="FEM")
ax[0].plot(fem_sol[:,0],fem_sol[:,1],color="blue",marker="o",markersize=5,label="FEM")
# Set labels for the axis
ax[0].set_xlabel(r"$x$")
ax[0].set_ylabel(r"$u(x)$")
ax[0].set_title("Verschiebungsverlauf")

# Set custom ticks for the x-axis with LaTeX symbols
tick_positions = [0, l/2, l]
tick_labels = [r"$0$", r"$\frac{\ell}{2}$", r"$\ell$"]
ax[0].set_xticks(tick_positions)
ax[0].set_xticklabels(tick_labels)
ax[0].set_yticks([0,F*l/(2*E*A)+3*n*l**2/(8*E),F*l/(E*A)+n*l**2 / (2*E)])
ax[0].set_yticklabels([r"$0$",r"$\frac{F\ell}{EA} + \frac{n \ell^2}{2E}$",r"$\frac{F\ell}{2EA} +\frac{3 n \ell^2}{8E}$"])

# Draw background grid
ax[0].grid(True)

# Add a second subplot (you can customize this as needed)
# For demonstration, let's just add a simple plot in the second subplot
fem_sol_N=np.array([F+3/4*n*A*l,F+1/4*n*A*l])
ax[1].plot(x, F+n*A*(l-x), color='red', linewidth=2.0)
ax[1].plot([0,l/2],[fem_sol_N[0],fem_sol_N[0]],color='blue',linewidth=2.0)
ax[1].plot([l/2,l],[fem_sol_N[1],fem_sol_N[1]],color='blue',linewidth=2.0)
ax[1].set_xlabel(r"$x$")
ax[1].set_ylabel(r"$N(x)$")
ax[1].set_title("Schnittkraftverlauf")
ax[1].grid(True)
ax[1].set_xticks(tick_positions)
ax[1].set_xticklabels(tick_labels)

ax[1].set_yticks([1,F+n*l/2*A,F+n*l])
ax[1].set_yticklabels([r"$F$",r"$F+\frac{nA\ell}{2}$",r"$F+nA\ell$"])

# Create a common legend below both subplots
fig.legend(loc='lower center', bbox_to_anchor=(0.5, -0.0), ncol=2)

# Adjust layout to make room for the legend
#plt.tight_layout(rect=[0, 0.1, 1, 1])  # Adjust the rect to make space for the legend

# Show the plot
plt.show()

# Add legend.
# ax[0].legend([r"analytische Lösung", r"FEM"])

from sympy import *
init_printing(use_unicode=True)

E,A,l,n = symbols('E A l n')
u_1, u_2, u_3 = symbols('u_1 u_2 u_3')
F, R = symbols('F R')

F_ext = Matrix([[F],[0],[R]])
F_v = (1/4)*n*l*A*Matrix([[1],[2],[1]])
F_g = F_ext+F_v
K=(2*E*A)/l*Matrix([[1,-1,0],[-1,2,-1],[0,-1,1]])
u=Matrix([[u_1],[u_2],[0]])

eq = F_ext+F_v - K*u
eq[0]

Kuu=K[0:2,0:2]
Kub=K[0:2,2]
Kbu=K[2,0:2]
Kbb=K[2,2]

F_g[0:2]

Kuu

rhs=Matrix(F_g[0:2])-Kub*u[2]
rhs

rhs=Matrix(F_g[0:2])-Kub*u[2]
display("rhs:",rhs)
uu = Kuu.LUsolve(rhs)
display("solution:",Kuu.LUsolve(rhs))

'rhs:'

'solution:'

Kbu*uu-Matrix([F_v[2]])

Ke=E*A/l*Matrix([[1,-1],[-1,1]])
Fve=(1/4)*n*l*A*Matrix([[1],[1]])
ue=Matrix([[uu[1]],[uu[0]]])
Ke*ue-Fve

eps=2*E*A/l * Matrix([[1, -1]])
N1=-eps*Matrix([[0],[uu[1]]])
N2=-eps*Matrix([[uu[1]],[uu[0]]])
display("N1",N1)
display("N2",N2)

'N1'

'N2'

uu

Stab plot

import numpy as np
import matplotlib.pyplot as plt

# Set up the plot configurations.
plt.rcParams['lines.linewidth'] = 2.0
plt.rcParams['lines.color'] = 'black'
plt.rcParams['legend.frameon'] = True
plt.rcParams['figure.figsize'] = (8, 6)
plt.rcParams['font.family'] = 'serif'
plt.rcParams['legend.fontsize'] = 15
plt.rcParams['font.size'] = 15
plt.rcParams['axes.spines.right'] = False
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.left'] = True
plt.rcParams['axes.spines.bottom'] = True
plt.rcParams['axes.axisbelow'] = True
plt.rcParams['grid.color'] = 'grey'
plt.rcParams['grid.alpha'] = 0.6
plt.rcParams['grid.linestyle'] = '--'

sizex = 7.5
sizey = 5

# Create a figure and axis.
fig, ax = plt.subplots(1, 2,figsize=(2*sizex,sizey))

F=1
E=1
A=1
l=20
n=1

x = np.linspace(0, l, 100)


# Plotting the functions.
ax[0].plot(x, (F/(E*A)+n*l/E)*x-(n/(2*E))*x**2,color='red',linewidth=2.0,label="analytische Lösung")

fem_sol=np.array([
 [0,0],
 [l/2,F*l/(2*E*A)+3*n*l**2/(8*E)],
 [l,F*l/(E*A)+n*l**2 / (2*E)]   
])

#ax[0].scatter(fem_sol[:,0],fem_sol[:,1],color="blue",marker="o",s=100,label="FEM")
ax[0].plot(fem_sol[:,0],fem_sol[:,1],color="blue",marker="o",markersize=5,label="FEM")
# Set labels for the axis
ax[0].set_xlabel(r"$x$")
ax[0].set_ylabel(r"$u(x)$")
ax[0].set_title("Verschiebungsverlauf")

# Set custom ticks for the x-axis with LaTeX symbols
tick_positions = [0, l/2, l]
tick_labels = [r"$0$", r"$\frac{\ell}{2}$", r"$\ell$"]
ax[0].set_xticks(tick_positions)
ax[0].set_xticklabels(tick_labels)
ax[0].set_yticks([0,F*l/(2*E*A)+3*n*l**2/(8*E),F*l/(E*A)+n*l**2 / (2*E)])
ax[0].set_yticklabels([r"$0$",r"$\frac{F\ell}{EA} + \frac{n \ell^2}{2E}$",r"$\frac{F\ell}{2EA} +\frac{3 n \ell^2}{8E}$"])

# Draw background grid
ax[0].grid(True)

# Add a second subplot (you can customize this as needed)
# For demonstration, let's just add a simple plot in the second subplot
fem_sol_N=np.array([F+3/4*n*A*l,F+1/4*n*A*l])
ax[1].plot(x, F+n*A*(l-x), color='red', linewidth=2.0)
ax[1].plot([0,l/2],[fem_sol_N[0],fem_sol_N[0]],color='blue',linewidth=2.0)
ax[1].plot([l/2,l],[fem_sol_N[1],fem_sol_N[1]],color='blue',linewidth=2.0)
ax[1].set_xlabel(r"$x$")
ax[1].set_ylabel(r"$N(x)$")
ax[1].set_title("Schnittkraftverlauf")
ax[1].grid(True)
ax[1].set_xticks(tick_positions)
ax[1].set_xticklabels(tick_labels)

ax[1].set_yticks([1,F+n*l/2*A,F+n*l])
ax[1].set_yticklabels([r"$F$",r"$F+\frac{nA\ell}{2}$",r"$F+nA\ell$"])

# Create a common legend below both subplots
fig.legend(loc='lower center', bbox_to_anchor=(0.5, -0.08), ncol=2)

# Adjust layout to make room for the legend
#plt.tight_layout(rect=[0, 0.1, 1, 1])  # Adjust the rect to make space for the legend

# Show the plot
plt.show()

# Add legend.
# ax[0].legend([r"analytische Lösung", r"FEM"])

Balken Plot

Shape Functions

import numpy as np
import matplotlib.pyplot as plt
import scienceplots

# Set up the plot configurations.
plt.rcParams['lines.linewidth'] = 2.0;
plt.rcParams['lines.color'] = 'black';
plt.rcParams['legend.frameon'] = True;
plt.rcParams['figure.figsize'] = (8, 6);
plt.rcParams['font.family'] = 'serif';
plt.rcParams['legend.fontsize'] = 15;
plt.rcParams['font.size'] = 15;
plt.rcParams['axes.spines.right'] = False;
plt.rcParams['axes.spines.top'] = False;
plt.rcParams['axes.spines.left'] = True;
plt.rcParams['axes.spines.bottom'] = True;
plt.rcParams['axes.axisbelow'] = True;
plt.rcParams['grid.color'] = 'grey';
plt.rcParams['grid.alpha'] = 0.6;
plt.rcParams['grid.linestyle'] = '--';

plt.style.use(['science', 'grid']);


fig,ax = plt.subplots(1,1);

# Define the x values
x = np.linspace(0, 1, 100);
# Define the shape functions - Hermite polynomials 4th order
N1 = 1 - 3*x**2 + 2*x**3;
N2 = x - 2*x**2 + x**3;
N3 = 3*x**2 - 2*x**3;
N4 = -x**2 + x**3;
# Plot the shape functions
ax.plot(x, N1, label='$N_1$');
ax.plot(x, N2, label='$N_2$');
ax.plot(x, N3, label='$N_3$');
ax.plot(x, N4, label='$N_4$');
# Add labels and legend
ax.set_xlabel(r'$\eta$');
ax.set_ylabel(r'$N(\eta)$');
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left');
# Show the plot
plt.show()

import numpy as np
import plotly.graph_objects as go

# Define the x values
x = np.linspace(0, 1, 100)

# Define the shape functions - Hermite polynomials 4th order
N1 = 1 - 3*x**2 + 2*x**3
N2 = x - 2*x**2 + x**3
N3 = 3*x**2 - 2*x**3
N4 = -x**2 + x**3

# Create the figure
fig = go.Figure()

# Add traces for each shape function
fig.add_trace(go.Scatter(x=x, y=N1, mode='lines', name='N_1'))
fig.add_trace(go.Scatter(x=x, y=N2, mode='lines', name='N_2'))
fig.add_trace(go.Scatter(x=x, y=N3, mode='lines', name='N_3'))
fig.add_trace(go.Scatter(x=x, y=N4, mode='lines', name='N_4'))

# Update layout with LaTeX math
fig.update_layout(
    title=r'$\text{Hermite Polynomials 4th Order}$',
    xaxis_title=r'$\eta$',
    yaxis_title=r'$N(\eta)$',
    legend=dict(x=1.05, y=1),
    font=dict(family="Serif", size=15),
    template='plotly_white',
    xaxis=dict(showgrid=True, gridcolor='grey', gridwidth=0.6, griddash='dash'),
    yaxis=dict(showgrid=True, gridcolor='grey', gridwidth=0.6, griddash='dash')
)

# Show the plot
fig.show()

import sympy as sp
from sympy import Matrix, latex

# Define the variables
xi,ell = sp.symbols('xi, ell')

# Define the Hermite polynomials
N1 = 1 - 3*xi**2 + 2*xi**3
N2 = (xi - 2*xi**2 + xi**3)*ell
N3 = 3*xi**2 - 2*xi**3
N4 = (-xi**2 + xi**3)*ell

N = Matrix([[N1], [N2], [N3], [N4]])
dNdxi = sp.diff(N, xi)
dNdxidxi = sp.diff(dNdxi, xi)

M=dNdxidxi*dNdxidxi.T
M

print(latex(M))

\left[\begin{matrix}\left(12 \xi - 6\right)^{2} & \ell \left(6 \xi - 4\right) \left(12 \xi - 6\right) & \left(6 - 12 \xi\right) \left(12 \xi - 6\right) & \ell \left(6 \xi - 2\right) \left(12 \xi - 6\right)\\\ell \left(6 \xi - 4\right) \left(12 \xi - 6\right) & \ell^{2} \left(6 \xi - 4\right)^{2} & \ell \left(6 - 12 \xi\right) \left(6 \xi - 4\right) & \ell^{2} \left(6 \xi - 4\right) \left(6 \xi - 2\right)\\\left(6 - 12 \xi\right) \left(12 \xi - 6\right) & \ell \left(6 - 12 \xi\right) \left(6 \xi - 4\right) & \left(6 - 12 \xi\right)^{2} & \ell \left(6 - 12 \xi\right) \left(6 \xi - 2\right)\\\ell \left(6 \xi - 2\right) \left(12 \xi - 6\right) & \ell^{2} \left(6 \xi - 4\right) \left(6 \xi - 2\right) & \ell \left(6 - 12 \xi\right) \left(6 \xi - 2\right) & \ell^{2} \left(6 \xi - 2\right)^{2}\end{matrix}\right]

print(latex(N))

\left[\begin{matrix}2 \xi^{3} - 3 \xi^{2} + 1\\\ell \left(\xi^{3} - 2 \xi^{2} + \xi\right)\\- 2 \xi^{3} + 3 \xi^{2}\\\ell \left(\xi^{3} - \xi^{2}\right)\end{matrix}\right]

K=sp.integrate(M, (xi, 0, 1))
K

print(latex(K))

\left[\begin{matrix}12 & 6 \ell & -12 & 6 \ell\\6 \ell & 4 \ell^{2} & - 6 \ell & 2 \ell^{2}\\-12 & - 6 \ell & 12 & - 6 \ell\\6 \ell & 2 \ell^{2} & - 6 \ell & 4 \ell^{2}\end{matrix}\right]

Nq = Matrix([[1-xi], [xi]])
q1,q2 = sp.symbols('q1, q2')
q = Matrix([[q1], [q2]])

Fq1 = N*(Nq.T*q)
Fq1

print(latex(Fq1))

\left[\begin{matrix}\left(q_{1} \left(1 - \xi\right) + q_{2} \xi\right) \left(2 \xi^{3} - 3 \xi^{2} + 1\right)\\\ell \left(q_{1} \left(1 - \xi\right) + q_{2} \xi\right) \left(\xi^{3} - 2 \xi^{2} + \xi\right)\\\left(- 2 \xi^{3} + 3 \xi^{2}\right) \left(q_{1} \left(1 - \xi\right) + q_{2} \xi\right)\\\ell \left(\xi^{3} - \xi^{2}\right) \left(q_{1} \left(1 - \xi\right) + q_{2} \xi\right)\end{matrix}\right]

Fqi = sp.integrate(Fq1, (xi, 0, 1))
Fqi

print(latex(Fqi))

\left[\begin{matrix}\frac{7 q_{1}}{20} + \frac{3 q_{2}}{20}\\\frac{\ell q_{1}}{20} + \frac{\ell q_{2}}{30}\\\frac{3 q_{1}}{20} + \frac{7 q_{2}}{20}\\- \frac{\ell q_{1}}{30} - \frac{\ell q_{2}}{20}\end{matrix}\right]

Balken unter konstantanter Streckenlast

from sympy import Rational
from sympy import lambdify
from sympy import Eq

x,EI,ell,q0,C_1,C_2,C_3,C_4 = sp.symbols('x,EI,ell,q0,C_1,C_2,C_3,C_4')

w = 1/EI*(Rational(1,24)*q0*x**4 +Rational(1,6)*C_1*x**3+Rational(1,2)*C_2*x**2+C_3*x+C_4)
wfun = lambdify(x, w)
w 

eq1 = Eq(w.subs(x,0),0)
eq2 = Eq(w.subs(x,ell),0)
eq3 = Eq(-w.diff(x,2).subs(x,0),0)
eq4 = Eq(-w.diff(x,2).subs(x,ell),0)
solution = sp.solve((eq1,eq2,eq3,eq4),(C_1,C_2,C_3,C_4))
solution

w.subs(solution)*EI

print(latex(w.subs(solution)*EI))

\frac{\ell^{3} q_{0} x}{24} - \frac{\ell q_{0} x^{3}}{12} + \frac{q_{0} x^{4}}{24}

print(latex(-w.subs(solution).diff(x,3)*EI))

q_{0} \left(\frac{\ell}{2} - x\right)

-w.subs(solution).diff(x,3)*EI

wp

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[33], line 1
----> 1 wp

NameError: name 'wp' is not defined

from sympy.plotting import plot, PlotGrid

wp = w.subs(solution).subs({"ell":1,"q0":1,"EI":1})
M = -wp.diff(x,2)
Q = -wp.diff(x,3)

p1= plot(-wp,(x,0,1),title='Durchbiegung $EI w$',show=False,xlabel=r'$\frac{x}{\ell}$',ylabel='$EI w$',line_color='black')
p2= plot(M,(x,0,1),title='Biegemoment $M$',show=False,xlabel=r'$\frac{x}{\ell}$',ylabel='$M$',line_color='red')
p3= plot(Q,(x,0,1),title='Querkraft $Q$',show=False,xlabel=r'$\frac{x}{\ell}$',ylabel='$Q$',line_color='blue')


PlotGrid(1,3,p1,p2,p3,size=(15,5))

<sympy.plotting.plotgrid.PlotGrid at 0x7f04740a5de0>

wc = w.subs(solution)
eq1 = Eq(wc.diff(x,1),0)
sol = sp.solve(eq1,x)
sol

[ell/2, ell/2 + sqrt(3)*ell/2, -sqrt(3)*ell/2 + ell/2]

wc.subs({"x":ell/2})

\[\displaystyle \frac{5 \ell^{4} q_{0}}{384 EI}\]